%% This function estimate the Variance-Covariance matrix of the vector theta of probability constraints

function V=vctp(Y,Z,D,B,A)

n=length(Y);

x=ceil(n*rand(n,B));
% Create bootstrap samples
yboot = Y(x);
dboot = D(x);
zboot = Z(x);
mu0=incosp(Y,Z,D,A);
mu=zeros(max(size(mu0)),B);
parfor i=1:B
    mu(:,i)=incosp(yboot(:,i),zboot(:,i),dboot(:,i),A);
    mu(:,i)=mu(:,i)-mu0;
end

Vb=cov(mu');

V=Vb*n;










